Linear Mixed Model

\[Y_{ij} = X_i\beta + Z_ja + \epsilon_{ij} \]

\[\epsilon \sim N(0, \sigma^2 \mathbf{I_n}) \] \[a \sim N(0, \sigma_a^2) \]

\[ y_i|a_{j} \sim𝑁(x_i \beta + a_j, \sigma^2) \]

\[ y_i \sim𝑁(x_i \beta, \sigma^2) \]

Standard Residual Plots

Standard Residuals

Raw \[\epsilon_i = Y_i - \hat{Y_i} \]

(internally) Studentized:

\[\epsilon_i = \frac {Y_i - \hat{Y_i}} {\hat{s}}\]

Pearson/Scaled:

\[ \epsilon_i = \frac {Y_i - \hat{Y_i}} {sd(Y_i)} \]

GLMM Formulation (binomial example)

Dependent variable: \(Y_{ij}\) (survived)

survival: \(p_i = Y_i/N\)

\[Y_{ij}|r_j \sim Binomial(N, \pi_{ij})\]

\[\eta_{ij} = \eta + \alpha_i + r_j\] (\(\eta\) = mean, \(\alpha_i\) = treatment, \(r_j\) = grouping variable)

\[\eta_{ij} = log(\frac{\pi_{ij}} {1 - \pi_{ij}} )\]

Standard Residuals from GLMMs Can Be Nonsensical

Raw/Scaled Residuals Can Lack Diagnostic Capabilities for GLMMs

• Visual patterns are difficult to interpret
• Over and under-dispersion are difficult to assess
• Goodness of fit tests are not valid

Background

For \(y_1,...,y_n\) responses:

\[y_i \sim \mathcal{P}(\mu_i, \phi)\]

CDF: \[F(y; \mu, \phi)\]

For a continuous \(F\), \(F(y_i; \mu_i, \phi)\) are uniformly distributed, and the quantile residuals are:

\[ r_{i} = \Phi^{-1} \left\{F(y_i; \hat{\mu_i}, \hat{\phi} )\right\}\]

Example: Survival Times of Patients

Survival: \(y_i \sim Exp(\mu_i)\)

\[ \text{log }\mu_i = \beta_0 + \beta_1 \text{ log } x_i \]

Quantile residuals:

\[ r_{i} = \Phi^{-1} \left\{ 1-exp(y_i/\hat{\mu_i}) \right\}\]

Discontinuous CDF, F

\[r_i = \Phi^{-1}(u_i)\] \[ u_i \sim Uniform(a_i, b_j]\] \[a_i = lim_{y \uparrow y_i} F(y; \hat{\mu_i}, \hat{\phi})\]

\[ b_i = F(y_i; \hat{\mu_i}, \hat{\phi}) \]

Binomial Example

\[y_i \sim binomial(n, p_i) \]

\[ \text{logit} (p_i) = \beta_0 + \beta_1x_i \] \[r_i = \Phi^{-1}(u_i)\]

\[ u_i \sim Uniform(a_i, b_j] \]

\[a_i = lim_{y \uparrow y_i} F(y; n, \hat{p_i}) =\displaystyle\sum_{i=0}^{\lfloor k-1 \rfloor} \begin{pmatrix} n \\ i \end{pmatrix} p_i(1-p)^{n-i}\]

\[b_i = F(y_i; n, \hat{p_i}) =\displaystyle\sum_{i=0}^{\lfloor k \rfloor} \begin{pmatrix} n \\ i \end{pmatrix} p_i(1-p)^{n-i}\]

Method implemented in ‘statmod’

This Method Has Been Around

The DHARMa Process

1. Model something using a generalized linear model with a chosen distribution and link function.

\[\mathbf{Y = X\beta+Za}\] \[ \mathbf{Y|a\sim G(\mu, R)} \] linear predictor: \(\mathbf{\eta = X\beta}\)

fitted value: \(E(Y) = \mu\)

Link function: \(\eta = g(\mu)\)

The DHARMa Process

For each observation in the data set:

2. Simulate new observations using fitted value as the mean and model parameters (e.g. shape, scale) for the other distributional parameters.

\[ Y_i \sim G(\hat{\mu_i}, \hat{\phi}) \]

3. Calculate the quantile for the cumulative distribution function

\[ r_{i} = F(y_i; \hat{\mu_i}, \hat{\phi} )\]

Expectations of the Quantile Residuals

\[ r_{ij} \sim Uniform(0,1) \]

\(r_{ij} = 0\): everything is larger

\(r_{ij} = 1\): everything is smaller

\(r_{ij} = 0.5\): right in the middle

DHARMa runs 250 simulations (per observation) by default, they recommend up to 1000

Poisson Example

\[ (\hat{Y_i} = \lambda = 5,; \quad Y_i = 7 )\]

Quantile Residuals for 1 Observation

Poisson GLMM Example

Create a data set with a random effect (“group”, 10 levels), a covariate (“Enviroinment1”) and a Poisson-distributed response variable (“observedResponse”):

testData = createData(sampleSize = 500)

Analyze correctly and incorrectly:

rightModel <- glmer(observedResponse ~ Environment1 + (1|group) , 
                     family = "poisson", data = testData)

wrongModel <- lmer(observedResponse ~ Environment1 + (1|group) , 
                     data = testData)

‘Standard’ Residuals

Quantile Residuals

Correctly Specified Model

Quantile Residuals

Incorrectly Specified Model

Distributions of Quantile Residuals

Quantile Function of the Standard Normal Distribution

Notes on DHARMa Implementation

• It depends on the simulation functions build into GLMM package (supported packages: lme4, glmmTMB, …)
• Commonly, only the last stochastic level (e.g. Poisson) is simulated, conditional on the fitted random effects - basically an omnibus test
• Residuals can be simulated for individuals predictors and tests can be conducted for individual factors
• for any kind of autocorrelation, conditional simulations are a must so that plots don’t indicate patterns or the tests fail due to autocorrelation

Conditional and Marginal Models

In Progress

(focus on lme4::glmer() and glmmTMB)

Other Functionality

(Implemented in DHARMa)

• Uniformity test
• Quantile location tesy
• Outlier test
• Tests for spatial and temporal autocorrelation
• Test for zero-inflation
• Over/Under dispersion tests

Best Practices

• Use quantile residuals because they are a better diagnostic feature for checking model fit
• Many of the tests are quite useful (but, as always, don’t rely on absolute threshold)
• Outlier test results should be treated with caution
• RTFM: read the DHARMa vignette; it is very helpful!
• ….More